Homework 14 – Small bacteria and ocean warming

Libraries

library(arm)
library(janitor)
library(car)
library(here)
library(tidyverse)
library(ggeffects)
library(MuMIn)
library(GGally)
library(gridExtra)
library(lme4)
library(lmerTest)
library(sjPlot)
library(RVAideMemoire)
library(DHARMa)
library(lubridate)
#install.packages("forecast")
library(forecast)
library(nlme)

Data

cell_raw <- read_csv(here("week_14", "data", "bacteria_size_temperature_edit.csv"))

LNA <- cell_raw %>%
  select(year,month,season,date,'day of year','temp 5 m E2','LNA ab uml','LNA bv','LNA B','%LNA BB') %>%
  drop_na()

LNA$date <- dmy(LNA$date)
LNA$season <- as.factor(LNA$season)
LNA$month <- as.factor(LNA$month)
LNA <- clean_names(LNA)

_‘LNA ab uml’ = LNA abundance in cells per mL, ‘LNA bv’ = LNA mean cell size in units of cubic microns, ‘LNA B’ = LNA total biomass in units of µg C per L, ‘%LNA BB’ = percent of total bacterial biomass in the LNA fraction.

Other important variables are year month date day of year, ‘temp 5 m E2’ = temperature at 5 meters depth._

1.

Begin by exploring seasonal patterns () in LNA abundance (LNA ab uml), cell size (LNA bv), biomass (LNA B), and percent biomass (%LNA BB), as well as temperature (temp 5m E2). I.e., make exploratory plots of these variables that focus on seasonality.

boxplot(lna_ab_uml ~ season, data = LNA, main = "LNA abundance", col = "lightgreen")
stripchart(lna_ab_uml ~ season, data = LNA,
           vertical = TRUE, method = "jitter",
           pch = 16, col = rgb(0, 0, 0, 0.5), add = TRUE)

boxplot(lna_bv ~ season, data = LNA, main = "LNA cell size", col = "lightblue")
stripchart(lna_bv ~ season, data = LNA,
           vertical = TRUE, method = "jitter",
           pch = 16, col = rgb(0, 0, 0, 0.5), add = TRUE)

boxplot(lna_b ~ season, data = LNA, main = "Biomass", col = "lightyellow")
stripchart(lna_b ~ season, data = LNA,
           vertical = TRUE, method = "jitter",
           pch = 16, col = rgb(0, 0, 0, 0.5), add = TRUE)

boxplot(percent_lna_bb ~ season, data = LNA, main = "% Biomass", col = "pink")
stripchart(percent_lna_bb ~ season, data = LNA,
           vertical = TRUE, method = "jitter",
           pch = 16, col = rgb(0, 0, 0, 0.5), add = TRUE)

boxplot(temp_5_m_e2 ~ season, data = LNA, main = "temperature at 5m depth", col = "gray")
stripchart(temp_5_m_e2 ~ season, data = LNA,
           vertical = TRUE, method = "jitter",
           pch = 16, col = rgb(0, 0, 0, 0.5), add = TRUE)

What have you learned so far?

LNA abundance might increase in the spring and summer

LNA cell size might increase in the winter and spring

Biomass might increase in the summer

The spring might be a period of disruption because percent biomass is all over the place

temperature at 5m deep is fairly consistent with seasonality

2.

Now make new exploratory plots of the same five variables, focusing on longer- term trends.

plot(LNA$date, LNA$lna_ab_uml, type = "l",xaxt = "n",
     main = "LNA abundance", xlab="",ylab="abund per ml")
monthly_ticks <- seq(from = as.Date("2002-04-18"), to = as.Date("2012-03-13"), by = "month")
axis(1, at = monthly_ticks, labels = format(monthly_ticks, "%b %Y"), las = 2,cex.axis=0.7)
lines(lowess(LNA$date, LNA$lna_ab_uml,f=0.1), col = "green",lwd = 2)

plot(LNA$date, LNA$lna_bv, type = "l",xaxt = "n",
     main = "LNA mean cell size", xlab="",ylab="cubic microns")
monthly_ticks <- seq(from = as.Date("2002-04-18"), to = as.Date("2012-03-13"), by = "month")
axis(1, at = monthly_ticks, labels = format(monthly_ticks, "%b %Y"), las = 2,cex.axis=0.7)
lines(lowess(LNA$date, LNA$lna_bv,f=0.1), col = "blue", lwd = 2)

plot(LNA$date, LNA$lna_b, type = "l",xaxt = "n",
     main = "Biomass", xlab="",ylab="cubic microns")
monthly_ticks <- seq(from = as.Date("2002-04-18"), to = as.Date("2012-03-13"), by = "month")
axis(1, at = monthly_ticks, labels = format(monthly_ticks, "%b %Y"), las = 2,cex.axis=0.7)
lines(lowess(LNA$date, LNA$lna_b,f=0.1), col = "purple", lwd = 2)

plot(LNA$date, LNA$percent_lna_bb, type = "l",xaxt = "n",
     main = "% Biomass", xlab="",ylab="Percent")
monthly_ticks <- seq(from = as.Date("2002-04-18"), to = as.Date("2012-03-13"), by = "month")
axis(1, at = monthly_ticks, labels = format(monthly_ticks, "%b %Y"), las = 2,cex.axis=0.7)
lines(lowess(LNA$date, LNA$percent_lna_bb,f=0.1), col = "red", lwd = 2)

plot(LNA$date, LNA$temp_5_m_e2, type = "l",xaxt = "n",
     main = "Temp", xlab="",ylab="Deg C")
monthly_ticks <- seq(from = as.Date("2002-04-18"), to = as.Date("2012-03-13"), by = "month")
axis(1, at = monthly_ticks, labels = format(monthly_ticks, "%b %Y"), las = 2,cex.axis=0.7)
lines(lowess(LNA$date, LNA$temp_5_m_e2,f=0.1), col = "black", lwd = 2)

How do you interpret these plots?

LNA abundance: peaks in August/September, 2008 to 2009 looks a little wonky (low throughout 2008, then very high around fall of 2009)

LNA cell size: peaks around April, “Crash” in January 2009 is wonky then does not recover (lower than other years)

Biomass: typically peaks in august/september, crash and large peak in Aug/Sep of 2009.

% Biomass: typically a large drop in April. April 2008 had substantially less of a dip, which carried on through at least 2012

Temp: fairly regular peaks and troughs where you’d expect

I am wondering if maybe there was, like, a bloom or something in 2008 and it takes a while to recover

3.

Let’s create statistical models to test for long-term trends, as well as relationships with temperature. Because the sampling is monthly, with only a couple gaps, we can treat the time series as one measured at discrete intervals. To keep track of the sample order you’ll need to construct a new column which numbers each sample in the order they were sampled. I.e., the column should look like 1, 2, 3, 4, etc., where the numbers correspond to the month in the time series. You can use the year and month columns to construct this.

LNA_sorted <- LNA[order(LNA$date), ]

LNA_sorted$sample_number <- seq_len(nrow(LNA_sorted))

Raw data: distributions

hist(log(LNA_sorted$lna_ab_uml))

LNA_sorted$ab_log <- log(LNA_sorted$lna_ab_uml)
hist(log(LNA_sorted$lna_bv))

LNA_sorted$bv_log<-log(LNA_sorted$lna_bv)
hist(sqrt(LNA_sorted$lna_b))

LNA_sorted$b_sq <-sqrt(LNA_sorted$lna_b)
hist(LNA_sorted$percent_lna_bb)

hist(LNA_sorted$temp_5_m_e2)

Use linear models to test for long term (linear) trends in the five variables: LNA abundance, LNA cell size, LNA biomass, LNA percent biomass, and temperature.

#abundance
lm.ab = gls(ab_log ~ sample_number, data = LNA_sorted,method = "ML")
LNA_sorted$resid.lmab.raw = resid(lm.ab)

ggAcf(LNA_sorted$resid.lmab.raw) + labs(title = 'abundance raw residuals, date as predictor')

#summary(lm.ab)
#cell size
lm.bv = gls(bv_log ~ sample_number, data = LNA_sorted,method = "ML")
LNA_sorted$resid.bv.raw = resid(lm.bv)
ggAcf(LNA_sorted$resid.bv.raw) + labs(title = 'cell size raw residuals, date predictor')

#summary(lm.bv)
#biomass 
lm.b = gls(b_sq ~ sample_number, data = LNA_sorted,method = "ML")
LNA_sorted$resid.b.raw = resid(lm.b)
ggAcf(LNA_sorted$resid.b.raw) + labs(title = 'biomass raw residuals, date predictor')

#summary(lm.b)
#percent biomass
lm.bb = gls(percent_lna_bb ~ sample_number, data = LNA_sorted,method = "ML")
LNA_sorted$resid.bb.raw = resid(lm.bb)
ggAcf(LNA_sorted$resid.bb.raw) + labs(title = 'percent biom raw residuals, date predictor')

#summary(lm.bb)
#temp
lm.t = gls(temp_5_m_e2 ~ sample_number, data = LNA_sorted,method = "ML")
LNA_sorted$resid.t.raw = resid(lm.t)
ggAcf(LNA_sorted$resid.t.raw) + labs(title = 'temp raw residuals, date predictor')

#summary(lm.t)

For each of these models, assess whether the residuals show evidence of temporal autocorrelation

It looks like they all do.

If autocorrelation is present, add a residual autocorrelation component to your model.

gl.ab = gls(ab_log ~ sample_number, data = LNA_sorted,correlation =corAR1(form =~ sample_number),method = "ML")
ggAcf(resid(gl.ab,type = 'normalized')) + labs(title = 'normalized residuals, abundance response')

#summary(gl.ab)
gl.bv = gls(bv_log ~ sample_number, data = LNA_sorted, correlation =corAR1(form =~ sample_number),method = "ML")
ggAcf(resid(gl.bv, type = 'normalized')) + labs(title = 'normalized residuals, cellsize response')

#summary(gl.bv)
gl.b = gls(b_sq ~ sample_number, data = LNA_sorted, correlation =corAR1(form =~ sample_number),method = "ML")
ggAcf(resid(gl.b, type = 'normalized')) + labs(title = 'normalized residuals, biomass response')

#summary(gl.b)
gl.bb = gls(percent_lna_bb ~ sample_number, data = LNA_sorted, correlation =corAR1(form =~ sample_number),method = "ML")
ggAcf(resid(gl.bb, type = 'normalized')) + labs(title = 'normalized residuals, %biomass response')

#summary(gl.bb)
gl.t = gls(temp_5_m_e2 ~ sample_number, data = LNA_sorted, correlation =corAR1(form =~ sample_number),method = "ML")
ggAcf(resid(gl.t, type = 'normalized')) + labs(title = 'normalized residuals, temp response')

#summary(gl.t)

Assess whether the model successfully accounted for autocorrelation.

…. not… uh, maybe? close for abundance, cell size, and biomass.less close for % biomass, not very close for temp. Was it supposed to? I feel like I’m not doing anything right anymore.

What have you learned from these models?

There is maybe some residual temporal autocorrelation, but there could be other factors explaining that autocorrelation too? I am also wondering if I should give up pretending I can learn and go work in one of Trump’s beautiful new tariff factories until I get made into glue like Boxer in Animal Farm I don’t know who cares wah wah baby crying

Does accounting for autocorrelation change the results? Yes, now just the temporal collection aspect isn’t enough. Which is probably more right than not for them all!

4.

Perform a similar set of analyses using temperature as the predictor, to test whether temperature can explain variation in the four LNA metrics, and whether autocorrelation needs to be accounted for when testing these relationships.

Abundance

gl.ab.t = gls(ab_log ~ temp_5_m_e2, data = LNA_sorted, method="ML")
ggAcf(resid(gl.ab.t, type = 'normalized')) + labs(title = 'normalized residuals, abundance response by temp')

#this one wins 



gl.ab.t.cor = gls(ab_log ~ temp_5_m_e2, data = LNA_sorted, correlation =corAR1(form =~ sample_number),method ="ML")
ggAcf(resid(gl.ab.t.cor, type = 'normalized')) + labs(title = 'normalized residuals, abundance response by temp AR1')

AICc(gl.ab.t)
## [1] 200.6047
AICc(gl.ab.t.cor)
## [1] 200.9486
summary(gl.ab.t)
## Generalized least squares fit by maximum likelihood
##   Model: ab_log ~ temp_5_m_e2 
##   Data: LNA_sorted 
##        AIC      BIC    logLik
##   200.3865 208.5951 -97.19327
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 11.168798 0.27744738 40.25555       0
## temp_5_m_e2  0.091266 0.01726342  5.28669       0
## 
##  Correlation: 
##             (Intr)
## temp_5_m_e2 -0.981
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -3.27447435 -0.56722140 -0.04596969  0.70474631  2.98658622 
## 
## Residual standard error: 0.5675843 
## Degrees of freedom: 114 total; 112 residual
summary(gl.ab.t.cor)
## Generalized least squares fit by maximum likelihood
##   Model: ab_log ~ temp_5_m_e2 
##   Data: LNA_sorted 
##        AIC      BIC    logLik
##   200.5817 211.5264 -96.29083
## 
## Correlation Structure: AR(1)
##  Formula: ~sample_number 
##  Parameter estimate(s):
##       Phi 
## 0.1303348 
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 11.213383 0.30373227 36.91864       0
## temp_5_m_e2  0.088255 0.01887816  4.67498       0
## 
##  Correlation: 
##             (Intr)
## temp_5_m_e2 -0.98 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -3.27843533 -0.55059328 -0.03466715  0.70402282  2.97993315 
## 
## Residual standard error: 0.5679109 
## Degrees of freedom: 114 total; 112 residual
#autocorrelation does not need to be accounted for

cell size

gl.bv.t = gls(bv_log ~ temp_5_m_e2, data = LNA_sorted, method="ML")
ggAcf(resid(gl.bv.t, type = 'normalized')) + labs(title = 'normalized residuals, cellsize response by temp')

#not good, need correlation


gl.bv.t.cor = gls(bv_log ~ temp_5_m_e2, data = LNA_sorted, correlation =corAR1(form =~ sample_number),method = "ML")
ggAcf(resid(gl.bv.t.cor, type = 'normalized')) + labs(title = 'normalized residuals, cellsize response by temp AR1')

#Better! Maybe need another level though 

gl.bv.tt = gls(bv_log ~ temp_5_m_e2, data = LNA_sorted, correlation =corARMA(form =~ sample_number, p = 3, q = 3),method = "ML")
ggAcf(resid(gl.bv.tt, type = 'normalized')) + labs(title = 'normalized residuals, cellsize response by temp AR2MA2')

#this looks pretty okay I think? 

AICc(gl.bv.t)
## [1] -159.2041
AICc(gl.bv.t.cor)
## [1] -172.6927
AICc(gl.bv.tt)
## [1] -179.015
summary(gl.bv.t)
## Generalized least squares fit by maximum likelihood
##   Model: bv_log ~ temp_5_m_e2 
##   Data: LNA_sorted 
##         AIC       BIC   logLik
##   -159.4223 -151.2137 82.71114
## 
## Coefficients:
##                  Value  Std.Error   t-value p-value
## (Intercept) -2.8531168 0.05725543 -49.83138  0.0000
## temp_5_m_e2 -0.0074707 0.00356257  -2.09700  0.0382
## 
##  Correlation: 
##             (Intr)
## temp_5_m_e2 -0.981
## 
## Standardized residuals:
##          Min           Q1          Med           Q3          Max 
## -3.479117512 -0.546169943  0.004079336  0.477672433  4.327886410 
## 
## Residual standard error: 0.1171295 
## Degrees of freedom: 114 total; 112 residual
summary(gl.bv.t.cor)
## Generalized least squares fit by maximum likelihood
##   Model: bv_log ~ temp_5_m_e2 
##   Data: LNA_sorted 
##         AIC       BIC   logLik
##   -173.0597 -162.1149 90.52985
## 
## Correlation Structure: AR(1)
##  Formula: ~sample_number 
##  Parameter estimate(s):
##       Phi 
## 0.3619635 
## 
## Coefficients:
##                  Value  Std.Error   t-value p-value
## (Intercept) -2.8511948 0.06994961 -40.76070  0.0000
## temp_5_m_e2 -0.0075344 0.00432522  -1.74198  0.0843
## 
##  Correlation: 
##             (Intr)
## temp_5_m_e2 -0.973
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -3.48530237 -0.55369128 -0.00259789  0.46830806  4.31490408 
## 
## Residual standard error: 0.1172485 
## Degrees of freedom: 114 total; 112 residual
summary(gl.bv.tt)
## Generalized least squares fit by maximum likelihood
##   Model: bv_log ~ temp_5_m_e2 
##   Data: LNA_sorted 
##         AIC     BIC  logLik
##   -180.7458 -156.12 99.3729
## 
## Correlation Structure: ARMA(3,3)
##  Formula: ~sample_number 
##  Parameter estimate(s):
##       Phi1       Phi2       Phi3     Theta1     Theta2     Theta3 
##  0.7631155 -0.2504339 -0.2230702 -0.5191559  0.4697360  0.3475245 
## 
## Coefficients:
##                  Value  Std.Error   t-value p-value
## (Intercept) -2.8135189 0.05991136 -46.96136  0.0000
## temp_5_m_e2 -0.0097503 0.00365270  -2.66935  0.0087
## 
##  Correlation: 
##             (Intr)
## temp_5_m_e2 -0.958
## 
## Standardized residuals:
##          Min           Q1          Med           Q3          Max 
## -3.454251489 -0.570368417 -0.003695439  0.451473959  4.121382481 
## 
## Residual standard error: 0.1213284 
## Degrees of freedom: 114 total; 112 residual

Biomass

gl.b.t = gls(b_sq ~ temp_5_m_e2, data = LNA_sorted, method="ML")
ggAcf(resid(gl.b.t, type = 'normalized')) + labs(title = 'normalized residuals, biomass response by temp')

#all good?


gl.b.t.cor = gls(b_sq ~ temp_5_m_e2, data = LNA_sorted, correlation =corAR1(form =~ sample_number),method = "ML")
ggAcf(resid(gl.b.t.cor, type = 'normalized')) + labs(title = 'normalized residuals, biomass response by temp AR1')

# not needed

AICc(gl.b.t)
## [1] 213.7405
AICc(gl.b.t.cor)
## [1] 213.1428
summary(gl.b.t)
## Generalized least squares fit by maximum likelihood
##   Model: b_sq ~ temp_5_m_e2 
##   Data: LNA_sorted 
##        AIC      BIC    logLik
##   213.5223 221.7309 -103.7611
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 0.8348747 0.29390136 2.840663  0.0053
## temp_5_m_e2 0.0874949 0.01828723 4.784482  0.0000
## 
##  Correlation: 
##             (Intr)
## temp_5_m_e2 -0.981
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.0523089 -0.7015614 -0.1651093  0.5635758  4.1830903 
## 
## Residual standard error: 0.6012448 
## Degrees of freedom: 114 total; 112 residual
summary(gl.b.t.cor)
## Generalized least squares fit by maximum likelihood
##   Model: b_sq ~ temp_5_m_e2 
##   Data: LNA_sorted 
##        AIC      BIC    logLik
##   212.7759 223.7206 -102.3879
## 
## Correlation Structure: AR(1)
##  Formula: ~sample_number 
##  Parameter estimate(s):
##       Phi 
## 0.1592262 
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept) 0.9222709 0.3275045 2.816055  0.0057
## temp_5_m_e2 0.0817733 0.0203484 4.018654  0.0001
## 
##  Correlation: 
##             (Intr)
## temp_5_m_e2 -0.979
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.0658441 -0.6785824 -0.1459732  0.5498167  4.1669360 
## 
## Residual standard error: 0.6016547 
## Degrees of freedom: 114 total; 112 residual

Percent Biomass

gl.bb.t = gls(percent_lna_bb ~ temp_5_m_e2, data = LNA_sorted, method="ML")
ggAcf(resid(gl.bb.t, type = 'normalized')) + labs(title = 'normalized residuals, %biomass response by temp')

#nope, autocorrelation


gl.bb.t.cor = gls(percent_lna_bb ~ temp_5_m_e2, data = LNA_sorted, correlation =corAR1(form =~ sample_number),method = "ML")
ggAcf(resid(gl.bb.t.cor, type = 'normalized')) + labs(title = 'normalized residuals, %biomass response by temp AR1')

#still a bit at lag 11

gl.bb.tt = gls(percent_lna_bb ~ temp_5_m_e2, data = LNA_sorted, correlation =corARMA(form =~ sample_number, p = 4, q = 4),method = "ML")
ggAcf(resid(gl.bb.tt, type = 'normalized')) + labs(title = 'normalized residuals, %biomass response 4AR4MA')

AICc(gl.bb.t)
## [1] -199.999
AICc(gl.bb.t.cor)
## [1] -204.631
AICc(gl.bb.tt)
## [1] -200.2229
summary(gl.bb.t)
## Generalized least squares fit by maximum likelihood
##   Model: percent_lna_bb ~ temp_5_m_e2 
##   Data: LNA_sorted 
##         AIC       BIC   logLik
##   -200.2172 -192.0086 103.1086
## 
## Coefficients:
##                  Value  Std.Error  t-value p-value
## (Intercept) 0.19678641 0.04787518 4.110406   1e-04
## temp_5_m_e2 0.01556887 0.00297891 5.226374   0e+00
## 
##  Correlation: 
##             (Intr)
## temp_5_m_e2 -0.981
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -3.24423144 -0.53826026  0.04110269  0.73018945  2.01331984 
## 
## Residual standard error: 0.09794002 
## Degrees of freedom: 114 total; 112 residual
summary(gl.bb.t.cor)
## Generalized least squares fit by maximum likelihood
##   Model: percent_lna_bb ~ temp_5_m_e2 
##   Data: LNA_sorted 
##         AIC       BIC  logLik
##   -204.9979 -194.0531 106.499
## 
## Correlation Structure: AR(1)
##  Formula: ~sample_number 
##  Parameter estimate(s):
##       Phi 
## 0.2482872 
## 
## Coefficients:
##                  Value  Std.Error  t-value p-value
## (Intercept) 0.21155723 0.05603314 3.775573   3e-04
## temp_5_m_e2 0.01456914 0.00347621 4.191096   1e-04
## 
##  Correlation: 
##             (Intr)
## temp_5_m_e2 -0.977
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -3.25743954 -0.50862611  0.06507804  0.72092525  2.00054770 
## 
## Residual standard error: 0.098116 
## Degrees of freedom: 114 total; 112 residual
summary(gl.bb.tt)
## Generalized least squares fit by maximum likelihood
##   Model: percent_lna_bb ~ temp_5_m_e2 
##   Data: LNA_sorted 
##         AIC      BIC   logLik
##   -202.8111 -172.713 112.4056
## 
## Correlation Structure: ARMA(4,4)
##  Formula: ~sample_number 
##  Parameter estimate(s):
##         Phi1         Phi2         Phi3         Phi4       Theta1       Theta2 
##  0.129090225  0.754155968 -0.020834376 -0.899247292 -0.004727259 -0.802165543 
##       Theta3       Theta4 
## -0.004483815  0.999749380 
## 
## Coefficients:
##                  Value  Std.Error  t-value p-value
## (Intercept) 0.20466903 0.07564865 2.705521  0.0079
## temp_5_m_e2 0.01506892 0.00475676 3.167898  0.0020
## 
##  Correlation: 
##             (Intr)
## temp_5_m_e2 -0.992
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -3.28066642 -0.53050205  0.04665563  0.72638741  2.01709447 
## 
## Residual standard error: 0.09728804 
## Degrees of freedom: 114 total; 112 residual

Make appropriate plots of the results and discuss the magnitude of the relationships.

#abundance 
ab<-plot(ggeffect(gl.ab.t, terms=c("temp_5_m_e2")))+ggtitle("LNA abundance")



#cell size
plot<-(ggeffect(gl.bv.tt, terms=c("temp_5_m_e2")))
plot$predicted <- exp(plot$predicted)
plot$conf.low <- exp(plot$conf.low)
plot$conf.high <- exp(plot$conf.high)
bv<-plot(plot)+ggtitle("LNA cell size")


#biomass 
b<-plot(ggeffect(gl.b.t, terms=c("temp_5_m_e2")))+ggtitle("Biomass")

#percent biomass 
bb<-plot(ggeffect(gl.bb.tt, terms=c("temp_5_m_e2")))+ggtitle("%Biomass")

grid.arrange(grobs=list(ab,bv,b,bb))

#summary(gl.ab.t)
#summary(gl.bv.tt)
#summary(gl.b.t)
#summary(gl.bb.tt)

Magnitudes from summaries (I know this isn’t very elegant and I could pull stuff out)

#Abundance
a<-exp(0.09127)-1
100*(a)
## [1] 9.556477
#Cell Size 
b<-exp(-0.0097503)-1
100*(b)
## [1] -0.970292
#Biomass 
c<-2*(0.0874949)
100*(c)
## [1] 17.49898
#%Biomass
100*(0.01506892)
## [1] 1.506892

Abundance: for every degree increase in C, abundance goes up by 9.6%

Cell size: for every degree increase in C, cell size goes DOWN by 0.97%

Biomass For every degree increase in C, biomass goes up by 17.5%

% Biomass For every degree increase in C, biomass goes up by 1.5%

Considering the whole set of analyses, what are your interpretations of the dynamics and drivers of LNA bacteria in this system?

I think when it gets hot LNA abundance and biomass go up, which just happens to be during the summer and is not so much dependent upon previous years/time periods– I think because cell size and percent biomass had some residual autocorrelation, they are more impacted by previous events. However, when things heat up cell size goes down just a bit, and biomass goes up by a bit. This makes sense looking at the exploratory time plots over the years.